!// Dmitry.Konovalov@jcu.edu.au Dmitry.A.Konovalov@gmail.com, 2013

module deriv_mod; 
use func_mod;
implicit none;
logical, save :: DERIV_TRACE = .FALSE.;
logical, save :: CACL_DRV_P9_TRACE = .TRUE.;

!from p24 of Bickley MathGaz 25 (1941)
! NOTE!!! Be careful, fortran mtrx is by-column
real(DP), dimension(5,9), save :: DRV_P9 = RESHAPE(& 
(/-109584._DP, 322560._DP, -564480._DP,  752640._DP, -705600._DP &
,   -5040._DP, -64224._DP,  141120._DP, -141120._DP,  117600._DP &
,     720._DP, -11520._DP,  -38304._DP,   80640._DP,  -50400._DP &
,    -240._DP,   2880._DP,  -20160._DP,  -18144._DP,   50400._DP &
,     144._DP,  -1536._DP,    8064._DP,  -32256._DP,        D0 &
,    -144._DP,   1440._DP,   -6720._DP,   20160._DP,  -50400._DP &
,     240._DP,  -2304._DP,   10080._DP,  -26880._DP,   50400._DP &
,    -720._DP,   6720._DP,  -28224._DP,   70560._DP, -117600._DP &
,    5040._DP, -46080._DP,  188160._DP, -451584._DP,  705600._DP/)&
, (/5,9/));

contains 


SUBROUTINE test_deriv_mod(); 
  character(LEN=64) :: name = "test_deriv_mod";
  call dbg_in(name);
  call numeric_deriv_test1();
  call dbg_out(name); 
end SUBROUTINE

SUBROUTINE numeric_deriv_test1();  
  integer :: nv;
  real(DP) :: minv, maxv, res, MAX_DRV_ERR = EPS18;
  real(DP), pointer :: px(:);      
  real(DP), dimension(200) :: py, py2, pdy, pf;      
  real(DP), dimension(3) :: coeff, dc;      
  real(DP), dimension(6) :: c6, dc6;      
  call dbg_str("-->numeric_deriv_test1()"); 

  nv = 9;  minv = D0; maxv = D1;
  px => mem_vec(nv);
  call grid_step(px,  minv, maxv);
  call dbg_assert_zero("9.px(1)=minv", px(1) - minv);
  call dbg_assert_zero("9.px(2)=minv", px(2) - (maxv+minv)/D8);
  call dbg_assert_zero("9.px(nv)=maxv", px(nv) - maxv);
  
  coeff = (/1, 0, 0/);  ! y=1
  call fun_vec_polynom(py, px, coeff);  
  call deriv_vec(pdy, px, py);
  dc = (/0, 0, 0/);  ! y=0
  call fun_vec_polynom(pf, px, dc);  
  call vec_assert_abs_diff_LE("pdy(1,0,0)", pdy, pf, MAX_DRV_ERR);

  coeff = (/D1, D1, D0/);  ! y=1+x 
  call fun_vec_polynom(py, px, coeff);  
  call deriv_vec(pdy, px, py);
  res = D1;
  call dbg_assert_abs_diff_LE("pdy(1)", pdy(1), res, MAX_DRV_ERR);
  call dbg_assert_abs_diff_LE("pdy(2)", pdy(2), res, MAX_DRV_ERR);
  call dbg_assert_abs_diff_LE("pdy(3)", pdy(3), res, MAX_DRV_ERR);
  call dbg_assert_abs_diff_LE("pdy(4)", pdy(4), res, MAX_DRV_ERR);
  call dbg_assert_abs_diff_LE("pdy(5)", pdy(5), res, MAX_DRV_ERR);
  call dbg_assert_abs_diff_LE("pdy(6)", pdy(6), res, MAX_DRV_ERR);
  call dbg_assert_abs_diff_LE("pdy(7)", pdy(7), res, MAX_DRV_ERR);
  call dbg_assert_abs_diff_LE("pdy(8)", pdy(8), res, MAX_DRV_ERR);
  call dbg_assert_abs_diff_LE("pdy(9)", pdy(9), res, MAX_DRV_ERR);

  coeff = (/D0, D0, D05/);  
  call fun_vec_polynom(py, px, coeff);  
  call deriv_vec(pdy, px, py);
  dc = (/D0, D1, D0/);  
  call fun_vec_polynom(pf, px, dc);  
  call vec_assert_abs_diff_LE("y=0.5*x^2", pdy, pf, MAX_DRV_ERR);

  coeff = (/D0, D0, -D05/);  
  call fun_vec_polynom(py, px, coeff);  
  call deriv_vec(pdy, px, py);
  dc = (/D0, -D1, D0/); 
  call fun_vec_polynom(pf, px, dc);  
  call vec_assert_abs_diff_LE("y=-0.5*x^2", pdy, pf, MAX_DRV_ERR);

  c6 = (/-D01, -D05, -D05, D0, D0, D0/);  
  call fun_vec_polynom(py, px, c6);  
  call deriv_vec(pdy, px, py);
  dc6 = (/-D05, -D1, D0, D0, D0, D0/);  
  call fun_vec_polynom(pf, px, dc6);  
  call vec_assert_abs_diff_LE("y=-0.1-0.5*x-0.5*x^2", pdy, pf, EPS14);

  c6 = (/-D02, D05, -D05, D2, -D1/D4, D0/);  
  call fun_vec_polynom(py, px, c6);  
  call deriv_vec(pdy, px, py);
  dc6 = (/D05, -D1, D6, -D1, D0, D0/);    
  call fun_vec_polynom(pf, px, dc6);  
  call vec_assert_abs_diff_LE("y=...x^4/4", pdy, pf, EPS14);

  call fun_vec_polynom(py, px, (/D0, D2/));  ! y=2x
  py2 = EXP(py);  ! y=exp(2x)
  call deriv_vec(pdy, px, py2);  ! y'=[exp(2x)]'=2exp(2x)  
  pf = py2 * 2;  
  call func_assert_abs_diff_LE("y=exp(2x)", px, pdy, pf, 2 * EPS5);

  ! NEW px
  nv = 101;  minv = -D1; maxv = D1;
  call mem_del_vec(px);
  px => mem_vec(nv);  
  call grid_step(px, minv, maxv); 
  call fun_vec_polynom(py, px, (/D0, D2/));  ! y=2x
  py2 = EXP(py);  ! y=exp(2x)
  call deriv_vec(pdy, px, py2);  ! y'=[exp(2x)]'=2exp(2x) 
  pf = py2 * 2;  
  call func_assert_abs_diff_LE("y=exp(2x)", px, pdy, pf, EPS10);

  c6 = (/-D02, D05, -D05, D2, -D1/D4, D0/);  
  call fun_vec_polynom(py, px, c6);  
  call deriv_vec(pdy, px, py);
  dc6 = (/D05, -D1, D6, -D1, D0, D0/);
  call fun_vec_polynom(pf, px, dc6);  
  call vec_assert_abs_diff_LE("y=...x^4/4", pdy, pf, EPS12);
  
  call mem_del_vec(px);
end SUBROUTINE


SUBROUTINE deriv_vec(res, px, py);    
  real(DP), intent(OUT) :: res(:);      
  real(DP), intent(IN) :: px(:), py(:);      
  call deriv_p9(res, px, py);
end SUBROUTINE

SUBROUTINE deriv_veca(res, px, funcA);    
  real(DP), intent(OUT) :: res(:,:);      
  real(DP), intent(IN) :: px(:), funcA(:,:);      
  integer i, mini, maxi, minf, maxf;
  !character(LEN=64) :: name = "deriv_veca";  
  !call dbg_in(name);

  ! make sure all dimensions are the same  
  call mtrx_i1_range(mini, maxi, res);
  call mtrx_i1_range(minf, maxf, funcA);
  call dbg_assert_true("1.mini==minf", mini == minf);
  call dbg_assert_true("1.maxi==maxf", maxi == maxf);

  call mtrx_i2_range(mini, maxi, res);
  call mtrx_i2_range(minf, maxf, funcA);
  call dbg_assert_true("2.mini==minf", mini == minf);
  call dbg_assert_true("2.maxi==maxf", maxi == maxf);

  do i = mini, maxi
    call deriv_vec(res(:,i), px, funcA(:,i));
  enddo
end SUBROUTINE

! 9-point derivative
! NOTE: this func fast but could be made faster if needed.
! Not tested for dimension(J:J2), where J != 1
SUBROUTINE deriv_p9(res, px, py);     
  real(DP), intent(OUT) :: res(:);      
  real(DP), intent(IN) :: px(:), py(:);      
  real(DP) stx, h2; ! step x
  integer xn, mini, maxi, i, k, fromi;
  integer yn, ymini, ymaxi;
  !call dbg_str("-->deriv_p9()", DERIV_TRACE);
  call vec_init_res(res, px);
  call vec_range(mini, maxi, px, xn);
  call vec_range(ymini, ymaxi, py, yn);
  call dbg_assert_true("xn >= 9", xn >= 9);
  call dbg_assert_true("ymini == mini", ymini == mini);
  call dbg_assert_true("maxi <= ymaxi", maxi <= ymaxi);
  call grid_step_assert_4(px);
  stx = px(mini+1) - px(mini);
  h2 = D05 / (20160D0 * stx);
 
  !call vec_dbg("px=", px);
  !call vec_dbg("py=", py);
  fromi = mini;
  do k = 1, 4 ! first 4
    i = fromi + k - 1;
    res(i) = h2 * calc_drv_p9(k, fromi, py);
    !call dbg_int("i=", i);
    !call dbg_dp("pdy(i)=", pdy(i));
  end do
  k = 5;
  do i = mini+4, maxi-4 ! middle
    fromi = i - k + 1;    ! i = fromi + k - 1;
    res(i) = h2 * calc_drv_p9(k, fromi, py);
    !call dbg_int("i=", i);
    !call dbg_dp("pdy(i)=", pdy(i));
  end do
  i = maxi-8;  
  do k = 6, 9 ! last 4
    i = fromi + k - 1;
    res(i) = h2 * calc_drv_p9(k, fromi, py);
    !call dbg_int("i=", (i+k-1));
    !call dbg_dp("pdy(i)=", pdy(i));
  end do
  
  return;
end SUBROUTINE

real(DP) function calc_drv_p9(k, fromi, py) result(res);     
  real(DP), intent(IN) :: py(:);      
  integer midi, k, fromi, mk;  
  logical :: trace;
  trace = CACL_DRV_P9_TRACE;
  !call dbg_int("k=", k, trace);
  !call dbg_int("fromi=", fromi, trace);
  midi = fromi + 4;   !call dbg_int("midi=", midi, trace);
  mk = 10 - k; !mirror k

  !call dbg_dp("DRV_P9(5, k)=", DRV_P9(5, k), trace);  
  res = DRV_P9(5, k) * py(midi); ! middle
  !call dbg_dp("res=", res, trace);  

  !call dbg_dp("  DRV_P9(4, k)=", DRV_P9(4, k), trace);  
  !call dbg_dp("DRV_P9(4, mk)=", DRV_P9(4, mk), trace);  
  res = res + DRV_P9(4,k)*py(midi-1) - DRV_P9(4,mk)*py(midi+1);
  !call dbg_dp("res=", res, trace);  

  !call dbg_dp("  DRV_P9(3, k)=", DRV_P9(3, k), trace);  
  !call dbg_dp("DRV_P9(3, mk)=", DRV_P9(3, mk), trace);  
  res = res + DRV_P9(3,k)*py(midi-2) - DRV_P9(3,mk)*py(midi+2);
  !call dbg_dp("res=", res, trace);  

  !call dbg_dp("  DRV_P9(2, k)=", DRV_P9(2, k), trace);  
  !call dbg_dp("DRV_P9(2, mk)=", DRV_P9(2, mk), trace);  
  res = res + DRV_P9(2,k)*py(midi-3) - DRV_P9(2,mk)*py(midi+3);
  !call dbg_dp("res=", res, trace);  

  res = res + DRV_P9(1,k)*py(midi-4) - DRV_P9(1,mk)*py(midi+4);
  !call dbg_dp("res=", res, trace);  
  return;
end function

end module 
